Trend Following System - Apple Stock

Calculate slow and fast exponential moving averages for AAPL stock using historical data, visualize results and calculate return and performance metrics. Our strategy purchases the stock as the fast exponential moving average crosses above the slow moving average. This system does not go short.

Alex Labuda (Suny New Paltz School of Business - Analytics Capstone)
9/9/2021

Retrieve Apple Stock Data

# silence warnings
getSymbols.warning4.0=FALSE

# stock data for Apple and the S&P 500
getSymbols("AAPL", src='yahoo', )
[1] "AAPL"
getSymbols("^GSPC", src='yahoo',)
[1] "^GSPC"
# Save apple data to a dataframe
df <- data.frame(Date=index(AAPL),coredata(AAPL))
# take most recent observations
df <- tail(df, 4000)
# reset index
rownames(df) = NULL

# SPX to a dataframe
SPX <- data.frame(Date=index(GSPC),coredata(GSPC))
# take most recent observations
SPX <- tail(SPX, 4000)
# reset index
rownames(SPX) = NULL
# copy the data
raw_data = df
# rename AAPL columns for simplicity
names(raw_data)[names(raw_data)=="AAPL.Open"] = "Open"
names(raw_data)[names(raw_data)=="AAPL.High"] = "High"
names(raw_data)[names(raw_data)=="AAPL.Low"] = "Low"
names(raw_data)[names(raw_data)=="AAPL.Close"] = "Close"
names(raw_data)[names(raw_data)=="AAPL.Volume"] = "Volume"
names(raw_data)[names(raw_data)=="AAPL.Adjusted"] = "Price"

# rename SPX columns for simplicity
names(SPX)[names(SPX)=="GSPC.Open"] = "Open"
names(SPX)[names(SPX)=="GSPC.High"] = "High"
names(SPX)[names(SPX)=="GSPC.Low"] = "Low"
names(SPX)[names(SPX)=="GSPC.Close"] = "Close"
names(SPX)[names(SPX)=="GSPC.Volume"] = "Volume"
names(SPX)[names(SPX)=="GSPC.Adjusted"] = "Price"

Sample Period

# retrieve max and min dates
date_max = max(raw_data$Date)
date_min = min(raw_data$Date)

print(paste("Our sample considers AAPL OHLC data from", date_min, "through", date_max))
[1] "Our sample considers AAPL OHLC data from 2007-01-03 through 2021-09-29"

Calculate Daily Return

Calculate a simple day-to-day return, from adjusted-close to adjusted-close:

# create a blank column
raw_data$return = NA
raw_data$log_return = NA

# loop through data and calculate daily returns
for(t in 2:nrow(raw_data)){
  raw_data$return[t] = (raw_data$Price[t] - raw_data$Price[t-1])/ raw_data$Price[t-1]
  raw_data$log_return[t] = log(raw_data$Price[t]/raw_data$Price[t-1])
}


# Do the same for SP500
SPX$return = NA
SPX$log_return = NA

for(t in 2:nrow(SPX)){
  SPX$return[t] = (SPX$Price[t] - SPX$Price[t-1])/ SPX$Price[t-1]
  SPX$log_return[t] = log(SPX$Price[t]/SPX$Price[t-1])
}

head(SPX, 4)
        Date    Open    High     Low   Close     Volume   Price
1 2007-01-03 1418.03 1429.42 1407.86 1416.60 3429160000 1416.60
2 2007-01-04 1416.60 1421.84 1408.43 1418.34 3004460000 1418.34
3 2007-01-05 1418.34 1418.34 1405.75 1409.71 2919400000 1409.71
4 2007-01-08 1409.26 1414.98 1403.97 1412.84 2763340000 1412.84
        return   log_return
1           NA           NA
2  0.001228286  0.001227532
3 -0.006084581 -0.006103168
4  0.002220318  0.002217857
head(raw_data,4)
        Date     Open     High      Low    Close     Volume    Price
1 2007-01-03 3.081786 3.092143 2.925000 2.992857 1238319600 2.569716
2 2007-01-04 3.001786 3.069643 2.993571 3.059286  847260400 2.626754
3 2007-01-05 3.063214 3.078571 3.014286 3.037500  834741600 2.608047
4 2007-01-08 3.070000 3.090357 3.045714 3.052500  797106800 2.620926
        return   log_return
1           NA           NA
2  0.022196227  0.021953476
3 -0.007121718 -0.007147198
4  0.004938178  0.004926025

Indicators and Studies

Simple Moving Averages Functions

Create a function to calculate two simple moving averages of differing lengths for use for a trend-following trading system:

Two_MovAvg_function = function(variable, slow_period = 100, fast_period = 20){
  # blank list of variables
  fast_ma = rep(NA, length(variable))
  slow_ma = rep(NA, length(variable))
  
  # loop to calculate average price for each record
  for (t in (slow_period+1):length(variable)){
    est_slow = mean(variable[(t-(slow_period+1)):(t-1)])
    est_fast = mean(variable[(t-(fast_period+1)):(t-1)])
    # update current average price
    fast_ma[t] = est_fast
    slow_ma[t] = est_slow
  }
  # create a dataframe for new vectors
  ma_data = data.frame(fast_ma = fast_ma,
                       slow_ma  = slow_ma)
  return(ma_data)
}

Exponential Moving Averages Function

Create a function to calculate two exponential moving averages of differing lengths for use for a trend-following trading system:

EA_Mov_Avg = function(variable, slow_lag = 100, fast_lag = 25){
  fast_ea = variable
  slow_ea = variable
  ema_diff = rep(NA, length(variable))
  
  # Loop to calculate each exponential average
  for (t in 2:length(variable)){ 
      est_slow = slow_ea[t-1] + (slow_ea[t] - slow_ea[t-1]) / ((slow_lag + 1) / 2)
      est_fast = fast_ea[t-1] + (fast_ea[t] - fast_ea[t-1]) / ((fast_lag + 1) / 2)
      
      ema_diff[t] = est_slow - est_fast
      
    slow_ea[t] = est_slow
    fast_ea[t] = est_fast
  }
  # dataframe for new vector
  ema_data = data.frame(fast_ema = fast_ea,
                        slow_ema = slow_ea,
                        ema_diff = ema_diff)
  return(ema_data)
}
# add these as new vectors to original data 

new_df = EA_Mov_Avg(raw_data$Close, slow_lag = 150, fast_lag = 30)

raw_data$slow_ma = new_df$slow_ema
raw_data$fast_ma = new_df$fast_ema
raw_data$ema_diff = new_df$ema_diff

Average True Range

True Range

First calculate the maximum of the aforementioned ranges:

max_range_functon = function(variable){
  max_range = rep(0, length(variable))
  
  for (t in 2:length(raw_data$Price)){
    max_rng = max(abs(raw_data$High[t] - raw_data$Low[t]), abs(raw_data$High[t] - raw_data$Close[t-1]), abs(raw_data$Close[t-1] - raw_data$Low[t]))
    
    max_range[t] = max_rng
  }
  max_range_data = data.frame(max_range = max_range)
  
  return(max_range_data)
}

# Test Function
# max_range_functon(raw_data$Close[1:10])

# Implement on full dataset
max_range_df = max_range_functon(raw_data$Close)

# Create new variable in data
raw_data$max_range = max_range_df$max_range
head(raw_data)
        Date     Open     High      Low    Close     Volume    Price
1 2007-01-03 3.081786 3.092143 2.925000 2.992857 1238319600 2.569716
2 2007-01-04 3.001786 3.069643 2.993571 3.059286  847260400 2.626754
3 2007-01-05 3.063214 3.078571 3.014286 3.037500  834741600 2.608047
4 2007-01-08 3.070000 3.090357 3.045714 3.052500  797106800 2.620926
5 2007-01-09 3.087500 3.320714 3.041071 3.306071 3349298400 2.838647
6 2007-01-10 3.383929 3.492857 3.337500 3.464286 2952880000 2.974494
        return   log_return  slow_ma  fast_ma     ema_diff max_range
1           NA           NA 2.992857 2.992857           NA  0.000000
2  0.022196227  0.021953476 2.993737 2.997143 -0.003405888  0.076786
3 -0.007121718 -0.007147198 2.994316 2.999746 -0.005429937  0.064285
4  0.004938178  0.004926025 2.995087 3.003150 -0.008062751  0.052857
5  0.083070258  0.079799840 2.999206 3.022693 -0.023487057  0.279643
6  0.047856250  0.046746410 3.005366 3.051183 -0.045816917  0.186786

Average True Range Function

Create function that uses the true range vector to calculate a moving average of the true range with a user-specified lag period:

ATR_function = function(variable, lag = 15){
  atr_list = rep(NA, length(variable))
  
  # Loop to create ATR
  for (t in (lag+2):length(variable)){
    est_atr = mean(variable[(t-(lag)):t])
    
    atr_list[t] = est_atr
  }
  
  atr_data = data.frame(atr = atr_list)
  return(atr_data)
}
# Add vector to our data
atr_data = ATR_function(raw_data$max_range)
raw_data$atr = atr_data$atr

Risk per lot

This is a multiple of the ATR for volatility-based position sizing:

risk_multiplier = 5
equity = 100000


raw_data$risk_per_lot = raw_data$atr * risk_multiplier
head(raw_data) 
        Date     Open     High      Low    Close     Volume    Price
1 2007-01-03 3.081786 3.092143 2.925000 2.992857 1238319600 2.569716
2 2007-01-04 3.001786 3.069643 2.993571 3.059286  847260400 2.626754
3 2007-01-05 3.063214 3.078571 3.014286 3.037500  834741600 2.608047
4 2007-01-08 3.070000 3.090357 3.045714 3.052500  797106800 2.620926
5 2007-01-09 3.087500 3.320714 3.041071 3.306071 3349298400 2.838647
6 2007-01-10 3.383929 3.492857 3.337500 3.464286 2952880000 2.974494
        return   log_return  slow_ma  fast_ma     ema_diff max_range
1           NA           NA 2.992857 2.992857           NA  0.000000
2  0.022196227  0.021953476 2.993737 2.997143 -0.003405888  0.076786
3 -0.007121718 -0.007147198 2.994316 2.999746 -0.005429937  0.064285
4  0.004938178  0.004926025 2.995087 3.003150 -0.008062751  0.052857
5  0.083070258  0.079799840 2.999206 3.022693 -0.023487057  0.279643
6  0.047856250  0.046746410 3.005366 3.051183 -0.045816917  0.186786
  atr risk_per_lot
1  NA           NA
2  NA           NA
3  NA           NA
4  NA           NA
5  NA           NA
6  NA           NA
# save moving average cross-over dates for chart annotations

h = raw_data[which((raw_data$fast_ma > raw_data$slow_ma) & (raw_data$fast_ma[-1] < raw_data$slow_ma[-1]))+1,]
l = raw_data[which((raw_data$fast_ma < raw_data$slow_ma) & (raw_data$fast_ma[-1] > raw_data$slow_ma[-1]))+1,]
h = na.exclude(h)
h
           Date     Open      High       Low     Close     Volume
276  2008-02-06  4.67250  4.711429  4.348929  4.357143 1573272400
429  2008-09-15  5.07250  5.274643  5.012857  5.012857  920634400
1479 2012-11-14 19.48214 19.551786 19.149286 19.174286  477170400
2168 2015-08-12 28.13250 28.855000 27.407499 28.809999  404870000
2997 2018-11-26 43.56000 43.737499 42.564999 43.654999  179994000
3336 2020-04-02 60.08500 61.287498 59.224998 61.232498  165934000
         Price      return  log_return   slow_ma  fast_ma   ema_diff
276   3.741115 -0.05689513 -0.05857779  5.524325  5.47893 0.04539427
429   4.304120 -0.05760740 -0.05933331  5.865623  5.82371 0.04191335
1479 16.609560 -0.01108855 -0.01115048 21.400512 21.37701 0.02350558
2168 26.453583  0.01541985  0.01530218 30.470711 30.39591 0.07480340
2997 42.441147  0.01352399  0.01343336 50.173922 49.85595 0.31796977
3336 60.567722  0.01668663  0.01654893 66.242964 66.13575 0.10721376
     max_range       atr risk_per_lot
276   0.362500 0.3355802     1.677901
429   0.306429 0.2209376     1.104688
1479  0.402500 0.6442408     3.221204
2168  1.447501 0.8901563     4.450782
2997  1.172500 1.7617186     8.808593
3336  2.062500 4.4723434    22.361717
l
           Date      Open      High       Low     Close     Volume
327  2008-04-21  5.793214  6.017857  5.777143  6.005714 1039152800
576  2009-04-16  4.256786  4.398214  4.242500  4.337500  593446000
1670 2013-08-20 18.203930 18.234644 17.886429 17.895357  358688400
2419 2016-08-10 27.177500 27.225000 26.940001 27.000000   96034000
3081 2019-03-29 47.457500 47.520000 47.134998 47.487499   94256000
3343 2020-04-14 70.000000 72.062500 69.512497 71.762497  194994800
         Price       return   log_return   slow_ma   fast_ma
327   5.156605  0.044212429  0.043262945  5.224567  5.236748
576   3.724249  0.032387180  0.031873771  3.858784  3.868324
1670 15.798877 -0.013137140 -0.013224196 16.457909 16.526846
2419 25.321693 -0.007444103 -0.007471948 25.482979 25.522523
3081 46.365067  0.006517443  0.006496297 45.244351 45.326908
3343 70.983398  0.050502962  0.049269060 66.259175 66.361778
         ema_diff max_range       atr risk_per_lot
327  -0.012180766  0.266428 0.1891517    0.9457584
576  -0.009539672  0.196785 0.1381247    0.6906237
1670 -0.068936203  0.348215 0.3600002    1.8000009
2419 -0.039543769  0.284999 0.4454686    2.2273431
3081 -0.082557051  0.385002 1.0609378    5.3046891
3343 -0.102603253  3.750000 3.0849998   15.4249991
# Create arrow annotations for trade entry signal

# Upper annotations
h_a <- list(
  x = h$Date,
  y = h$slow_ma,
  text = "Sell",
  xref = "x",
  yref = "y",
  showarrow = TRUE,
  arrowcolor = "red",
  arrowhead = 5,
  arrowsize = 0.8,
  arrowwidth = 1.5,
  opacity = 0.75,
  align = "left",
  ax = 5,
  ay = -55
)

# Lower annotations

l_a <- list(
  x = l$Date,
  y = l$slow_ma,
  text = "Buy",
  xref = "x",
  yref = "y",
  showarrow = TRUE,
  arrowcolor = "green",
  arrowhead = 5,
  arrowsize = 0.8,
  arrowwidth = 1.5,
  opacity = 0.75,
  align = "right",
  ax = -5,
  ay = 55
)


# figure labels
f <- list(
  family = "Courier New, monospace",
  size = 18,
  color = "#7f7f7f ")
x <- list(
  title = "x Axis",
  titlefont = f)
y <- list(
  title = "y Axis",
  titlefont = f)

Candlestick Chart

Visualizing the Trend-following System

# colors column for increasing and decreasing
for (i in 2:length(raw_data[,1])) {
  if (raw_data$ema_diff[i] > 0) {
      raw_data$MACD_direction[i] = 'Increasing'
  } else {
      raw_data$MACD_direction[i] = 'Decreasing'
  }
}

# color of bars for chart
i <- list(line = list(color = '#17BECF')) # 'i' for increasing
d <- list(line = list(color = '#b22222')) # 'd' for decreasing


# Plot main data
fig <- raw_data %>% plot_ly(x = ~Date, type="ohlc",
          open = ~Open, close = ~Close,
          high = ~High, low = ~Low, name = "AAPL",
          increasing = i, decreasing = d)

# Add Fast and Slow moving average lines
fig <- fig %>% add_lines(x = ~Date, y = ~slow_ma, name = "Slow",
            line = list(color = '#ccc', width = 0.5),
            legendgroup = "Bands", inherit = F,
            showlegend = TRUE, hoverinfo = "none") 
fig <- fig %>% add_lines(x = ~Date, y = ~fast_ma, name = "Fast",
            line = list(color = '#E377C2', width = 0.5),
            hoverinfo = "none", inherit = F)
# Add y-axis title
fig <- fig %>% layout(yaxis = list(title = "Price"))


# plot MACD line subplot
fig2 <- raw_data 
fig2 <- fig2 %>% plot_ly(x=~Date, y=~ema_diff, type='bar', name = "MACD",
                         color = ~MACD_direction, colors = c("red", "green"), alpha = 0.8)
fig2 <- fig2 %>% layout(yaxis = list(title = "MACD"))

# use subplot to add plots to the same figure
fig = subplot(fig, fig2, heights = c(0.8,0.15), nrows = 2,
              shareX = TRUE, titleY = TRUE)

# Add arrow annotations
fig <- fig %>% layout(annotations = h_a)
fig <- fig %>% layout(annotations = l_a)

# Add title
fig <- fig %>% layout(title = "AAPL Stock - Moving Average Crossover",
                       xaxis = list(rangeslider = list(visible = F)))

fig

Generate signals using MA Indicator

Now we use our indicator to generate trading signals and store them in a vector:

raw_data$signal = NA

# Loop through data to create signal column
for(t in 100:nrow(raw_data)){
  
  current_slow = raw_data$slow_ma[t]
  current_fast = raw_data$fast_ma[t]
  
  prev_slow = raw_data$slow_ma[t-1]
  prev_fast = raw_data$fast_ma[t-1]
  
  # use 'if' statement to translate crossover into signals
  if(current_fast > current_slow & prev_fast < prev_slow) {signal = "B"} else 
    if(current_fast < current_slow & prev_fast > prev_slow) {signal = "S"} else
      {signal = "H"}
  
  raw_data$signal[t] = signal
}

Holding Status and Investment Return

After generating the signals, we must calculate the return from our strategy:

# used to calculate fill price penalty / skid
# we will penalize by 40% from the open - high for buy orders
# and 40% from the open - low for sell orders

buy_weight = .40
sell_weight = -.40
# Generate a holding status column  
# and investment return column

raw_data$holding = 0      # first record, you are not holding
raw_data$Inv_return = NA  # obviously the return will also be 0

# we can begin the loop
for (t in 100:nrow(raw_data)){
  
  if(t==1){prev_holding=0} 
  else {prev_holding=raw_data$holding[t-1]}
  
  # look at the signal to decide the change of holding status
  signal = raw_data$signal[t]
  if(signal=="B"){raw_data$holding[t]=1} else 
    if (signal=="S") {raw_data$holding[t]=0} else 
      if (signal=="H"){raw_data$holding[t]=prev_holding}
  
  # Now calculate investment return
  Open = raw_data$Open[t]
  Close = raw_data$Close[t]
  High = raw_data$High[t]
  Low = raw_data$Low[t]
  
  if(prev_holding==0 & signal=="B"){investment_return=Close/((High-Open)*buy_weight+Open) - 1} else 
    if(prev_holding==0 & signal=="H"){investment_return=0} else 
      #if(prev_holding==1 & signal=="S"){investment_return=Close[t-1]/(Open*(1+sell_weight)) - 1} else
      #if(prev_holding==0 & signal=="S"){investment_return=(-1)*(Close/(Open*(1+sell_weight)) - 1)} else
        if(prev_holding==1 & signal=="B"){investment_return=Close/raw_data$Close[t-1]-1} else
          if(prev_holding==1 & signal=="H"){investment_return=Close/raw_data$Close[t-1]-1} else
            {investment_return = Open/raw_data$Close[t-1]-1}
 
   # save the investment return to the dataset
  raw_data$Inv_return[t] = investment_return}

Cumulative Return

From the individual daily returns, we must create a cash and number of stocks vector to hold each days current figures:

# starting investment $100,000
initial_wealth = 100000
# create vectors to store updates
raw_data$cash = rep(initial_wealth, length(raw_data$Close))
raw_data$n_stock = rep(0, length(raw_data$Close))

# loop through signals and simulate the trades
for (t in 1:length(raw_data$signal)){
  
  if(t==1){prev_cash = initial_wealth; n_stock = 0
    
    } else
      
      {prev_cash = raw_data$cash[t-1]; 
      n_stock = raw_data$n_stock[t-1]}
  
  signal = raw_data$signal[t]
  Open = raw_data$Open[t]
  Close = raw_data$Close[t]
  High = raw_data$High[t]
  Low = raw_data$Low[t]
  
# is.na to deal with rows at  beginning with no signal
  
  if(is.na(signal)==TRUE){
    prev_cash = raw_data$cash[t-1]
    n_stock = raw_data$n_stock[t-1] 
    
    } else
      
      if (signal=="B"){
        buy_skid = (High-Open) * buy_weight + Open
        trans_prc = buy_skid
        n_change = floor(prev_cash/trans_prc)
        stock_n = n_stock + n_change
        cash = prev_cash - n_change * trans_prc
            
        raw_data$cash[t] = cash
        raw_data$n_stock[t] = stock_n
        
        } else
          
          if (signal =="S"){
            sell_skid = (Open-Low) * sell_weight + Open
            trans_prc = sell_skid
            n_change = n_stock
            stock_n = n_change - n_stock
            cash = prev_cash + n_change * trans_prc
        
            raw_data$cash[t] = cash
            raw_data$n_stock[t] = stock_n
            
            }else
              
              if (signal == "H") {
                raw_data$cash[t] = raw_data$cash[t-1]
                raw_data$n_stock[t] = raw_data$n_stock[t-1]
                
                }
}

Final Return

Calculate the cash plus the value of the stock on hand for the last day in the sample:

print(raw_data[max(index(raw_data)),])
           Date   Open   High    Low  Close   Volume  Price
3712 2021-09-29 142.47 144.45 142.03 142.83 74602000 142.83
          return  log_return  slow_ma  fast_ma  ema_diff max_range
3712 0.006482968 0.006462044 137.9734 147.2232 -9.249815  2.539993
          atr risk_per_lot MACD_direction signal holding  Inv_return
3712 3.065625     15.32813     Decreasing      H       1 0.006482968
         cash n_stock
3712 60.70275   16587
final_return = raw_data$n_stock[max(index(raw_data))] * raw_data$Close[max(index(raw_data))] + raw_data$cash[max(index(raw_data))-1]

print(paste("Initial Investment:", sprintf("$ %3.2f", initial_wealth)))
[1] "Initial Investment: $ 100000.00"
print(paste("Final Return:",sprintf("$ %3.2f" , final_return)))
[1] "Final Return: $ 2369181.95"

Visualize Days with Largest Loss and Gain

max_gain=raw_data[which.max(raw_data$Inv_return),]
max_loss = raw_data[which.min(raw_data$Inv_return),]

# max(raw_data$Inv_return, na.rm = TRUE)

# Create subset
max_data = raw_data[3310:3330,]
min_data = raw_data[3310:3330,]
max_loss
           Date    Open  High Low   Close    Volume    Price
3323 2020-03-16 60.4875 64.77  60 60.5525 322423600 59.89511
         return log_return  slow_ma  fast_ma  ema_diff max_range
3323 -0.1286471 -0.1377082 67.15978 72.53165 -5.371863    9.4925
          atr risk_per_lot MACD_direction signal holding Inv_return
3323 5.212187     26.06094     Decreasing      H       1 -0.1286469
         cash n_stock
3323 28.24028   19665
max_gain
           Date    Open  High     Low   Close    Volume    Price
3322 2020-03-13 66.2225 69.98 63.2375 69.4925 370732000 68.73806
        return log_return  slow_ma  fast_ma  ema_diff max_range
3322 0.1198085  0.1131577 67.24847 73.35779 -6.109322  7.922504
          atr risk_per_lot MACD_direction signal holding Inv_return
3322 4.774375     23.87188     Decreasing      H       1  0.1198083
         cash n_stock
3322 28.24028   19665
# Create arrow annotations

# Upper annotations
h_a <- list(
  x = max_gain$Date,
  y = max_gain$Price,
  text = "Largest Gain",
  xref = "x",
  yref = "y",
  showarrow = TRUE,
  arrowcolor = "green",
  arrowhead = 5,
  arrowsize = 0.8,
  arrowwidth = 2.5,
  opacity = 0.75,
  align = "right",
  ax = 5,
  ay = -55
)

# Lower annotations

l_a <- list(
  x = max_loss$Date,
  y = max_loss$Price,
  text = "Largest Loss",
  xref = "x",
  yref = "y",
  showarrow = TRUE,
  arrowcolor = "green",
  arrowhead = 5,
  arrowsize = 0.8,
  arrowwidth = 2.5,
  opacity = 0.75,
  align = "right",
  ax = -5,
  ay = 55
)


# figure labels
f <- list(
  family = "Courier New, monospace",
  size = 18,
  color = "#7f7f7f ")
x <- list(
  title = "x Axis",
  titlefont = f)
y <- list(
  title = "y Axis",
  titlefont = f)

Largest Gain

fig <- max_data %>% plot_ly(x = ~Date, type="ohlc",
          open = ~Open, close = ~Close,
          high = ~High, low = ~Low, name = "AAPL",
          increasing = i, decreasing = d)

# color of bars for chart
i <- list(line = list(color = '#17BECF')) # 'i' for increasing
d <- list(line = list(color = '#b22222')) # 'd' for decreasing

# Add Fast and Slow moving average lines
fig <- fig %>% add_lines(x = ~Date, y = ~slow_ma, name = "Slow",
            line = list(color = '#ccc', width = 0.5),
            legendgroup = "Bands", inherit = F,
            showlegend = TRUE, hoverinfo = "none") 
fig <- fig %>% add_lines(x = ~Date, y = ~fast_ma, name = "Fast",
            line = list(color = '#E377C2', width = 0.5),
            hoverinfo = "none", inherit = F)

# Add y-axis title
fig <- fig %>% layout(yaxis = list(title = "Price"))

# Add arrow annotations
fig <- fig %>% layout(annotations = h_a)

# Add title
fig <- fig %>% layout(title = "AAPL Stock - Strategy Largest 1-Day Gain",
                       xaxis = list(rangeslider = list(visible = F)))

fig

Largest Loss

fig <- min_data %>% plot_ly(x = ~Date, type="ohlc",
          open = ~Open, close = ~Close,
          high = ~High, low = ~Low, name = "AAPL",
          increasing = i, decreasing = d)

# color of bars for chart
i <- list(line = list(color = '#17BECF')) # 'i' for increasing
d <- list(line = list(color = '#b22222')) # 'd' for decreasing

# Add Fast and Slow moving average lines
fig <- fig %>% add_lines(x = ~Date, y = ~slow_ma, name = "Slow",
            line = list(color = '#ccc', width = 0.5),
            legendgroup = "Bands", inherit = F,
            showlegend = TRUE, hoverinfo = "none") 
fig <- fig %>% add_lines(x = ~Date, y = ~fast_ma, name = "Fast",
            line = list(color = '#E377C2', width = 0.5),
            hoverinfo = "none", inherit = F)

# Add y-axis title
fig <- fig %>% layout(yaxis = list(title = "Price"))

# Add arrow annotations
fig <- fig %>% layout(annotations = l_a)

# Add title
fig <- fig %>% layout(title = "AAPL Stock - Strategy Largest 1-day Loss",
                       xaxis = list(rangeslider = list(visible = F)))

fig

Max Adverse & Favorable Single-Day Returns

max_daily_return = round((max(raw_data$Inv_return, na.rm = TRUE)*100),4)
max_daily_loss =  round((min(raw_data$Inv_return, na.rm = TRUE)*100), 4)
avg_daily_return = round((mean(raw_data$Inv_return, na.rm = TRUE)*100),4)

bh_max_return = round((max(raw_data$return, na.rm = TRUE)*100),4)
bh_max_loss =  round((min(raw_data$return, na.rm = TRUE)*100), 4)
bh_avg_return = round((mean(raw_data$return, na.rm = TRUE)*100),4)

print(paste("Largest 1-day return with trend system:", max_daily_return,"%"))
[1] "Largest 1-day return with trend system: 11.9808 %"
print(paste("Largest 1-day loss with trend system:", max_daily_loss,"%"))
[1] "Largest 1-day loss with trend system: -12.8647 %"
print(paste("Average 1-day with trend system return:", avg_daily_return,"%"))
[1] "Average 1-day with trend system return: 0.1004 %"
print(paste("Largest 1-day return buy and hold strategy:", bh_max_return,"%"))
[1] "Largest 1-day return buy and hold strategy: 13.905 %"
print(paste("Largest 1-day loss with buy and hold strategy:", bh_max_loss,"%"))
[1] "Largest 1-day loss with buy and hold strategy: -17.9195 %"
print(paste("Average 1-day with buy and hold strategy:", bh_avg_return,"%"))
[1] "Average 1-day with buy and hold strategy: 0.129 %"

Performance Evaluation

Distribution of Returns

Create a histogram of returns for a simple buy and hold strategy, as well as for our trend-following system:

p <- raw_data %>%
  ggplot(aes(x=Inv_return)) +
    geom_histogram(fill="#69b3a2", color="#e9ecef", alpha=0.7, binwidth = .005) +
  coord_cartesian(xlim = c(-0.07, 0.07), ylim = c(0, 500)) +
  labs(x = "Investment Returns", y = "Frequency of Returns",
       title = "EMA Crossover - Histogram of Investment Returns",
       caption = "Trend following system - Slow lag: 100 periods, Fast lag: 20 periods")+
  theme_classic()
p
bh = raw_data %>%
  ggplot(aes(x=return)) +
    geom_histogram(fill="#69b3a2", color="#e9ecef", alpha=0.7, binwidth = .005) +
  coord_cartesian(xlim = c(-0.07, 0.07), ylim = c(0, 500)) +
  labs(x = "Investment Returns", y = "Frequency of Returns",
       title = "Buy / Hold - Histogram of Investment Returns",
       caption = "This histogram considers a buy and hold strategy")+
  theme_classic()
bh

Mean and Standard Deviation

Retrieve the average daily return and volatility of our trend-following system, as well as a simple buy and hold strategy:

average_return = mean(raw_data$Inv_return, na.rm = TRUE) # average return
volatility = sd(raw_data$Inv_return, na.rm = TRUE) # volatility

print(paste("Average return:", round(average_return,6)))
[1] "Average return: 0.001004"
print(paste("Volatility:", round(volatility,6)))
[1] "Volatility: 0.015134"
bh_average_return = mean(raw_data$return, na.rm = TRUE) # average return
bh_volatility = sd(raw_data$return, na.rm = TRUE) # volatility

print(paste("Average return:", round(bh_average_return,6)))
[1] "Average return: 0.00129"
print(paste("Volatility:", round(bh_volatility,6)))
[1] "Volatility: 0.020302"
print(paste("Buy and hold strategy experiences much higher volatility for a small amount of additional return"))
[1] "Buy and hold strategy experiences much higher volatility for a small amount of additional return"

Sharpe Ratio

The Sharpe ratio was developed by Nobel laureate William F. Sharpe and is used to help investors understand the return of an investment compared to its risk.12 The ratio is the average return earned in excess of the risk-free rate per unit of volatility or total risk. Volatility is a measure of the price fluctuations of an asset or portfolio.

rf = .00001 # risk-free rate
print(paste("Trend following: Sharpe ratio of daily returns:",round((average_return-rf)/volatility,5)))
[1] "Trend following: Sharpe ratio of daily returns: 0.0657"
print(paste("Buy & hold: Sharpe ratio of daily returns:",round((bh_average_return-rf)/bh_volatility,5)))
[1] "Buy & hold: Sharpe ratio of daily returns: 0.06304"
print(paste("Our Trend-following system achieves slightly higher Sharpe ratio"))
[1] "Our Trend-following system achieves slightly higher Sharpe ratio"

Value at Risk

VAR measures the cut-off return that your financial asset will fall below with certain probability.

Value at risk (VaR) is a statistic that quantifies the extent of possible financial losses within a firm, portfolio, or position over a specific time frame. This metric is most commonly used by investmentand commercial banks to determine the extent and probabilities of potential losses in their institutional portfolios.

var = quantile(raw_data$Inv_return, 0.05, na.rm = TRUE)*100
bh_var = quantile(raw_data$return, 0.05, na.rm = TRUE)*100

print(paste("Trend-following VAR: we can expect 95% of the time returns will be greater than:", round(var,4)))
[1] "Trend-following VAR: we can expect 95% of the time returns will be greater than: -2.361"
print(paste("Buy / Hold model VAR: we can expect 95% of the time returns will be greater than:", round(bh_var,4)))
[1] "Buy / Hold model VAR: we can expect 95% of the time returns will be greater than: -2.9629"

Summary of Investment & Market Return

In order to break down risk into two components, we need to collect what part of the return was attributed to the market and what part to our trading system.

Fit CAPM Model

y = raw_data$Inv_return - rf
x = SPX$return - rf
# fit regression

capm = lm(y~x)
summary(capm)

Call:
lm(formula = y ~ x)

Residuals:
      Min        1Q    Median        3Q       Max 
-0.084232 -0.006302 -0.000674  0.006462  0.099507 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) 0.0007829  0.0002192   3.572  0.00036 ***
x           0.5729293  0.0168314  34.039  < 2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.01317 on 3611 degrees of freedom
  (99 observations deleted due to missingness)
Multiple R-squared:  0.2429,    Adjusted R-squared:  0.2427 
F-statistic:  1159 on 1 and 3611 DF,  p-value: < 2.2e-16
print(paste("Trend-following system Beta:", round(capm$coefficients[2],4)))
[1] "Trend-following system Beta: 0.5729"

The larger the \(BETA\) is, the larger the systematic risk is.

We use this to measure the aggressiveness of the stock or the trading system.

Our trend-following system is mostly conservative.

Jensen’s Alpha (Abnormal Return)

tf_beta = capm$coefficients[2] # beta
jensens_alpha = capm$coefficients[1] # Jensen's Alpha

print(paste("Jensen's Alpha:", round(jensens_alpha,4)))
[1] "Jensen's Alpha: 8e-04"

Treynor’s Ratio

print(paste("Treynor's Ratio Version 1:", round(jensens_alpha/tf_beta, 4)))
[1] "Treynor's Ratio Version 1: 0.0014"
print(paste("Treynor's Ratio Version 2:", round(((mean(raw_data$Inv_return, na.rm = TRUE))/tf_beta),4)))
[1] "Treynor's Ratio Version 2: 0.0018"